Inhomogeneous Fermi mixtures at Unitarity: 
Bogoliubov-de Gennes vs. Landau-Ginzburg 
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We present an inhomogeneous theory for the low-temperature properties of a resonantly inter- 
acting Fermi mixture in a trap that goes beyond the local-density approximation. We compare 
the Bogoliubov-de Gennes and a Landau-Ginzburg approach and conclude that the latter is more 
appropriate when dealing with a first-order phase transition. Our approach incorporates the state- 
of-the-art knowledge on the homogeneous mixture with a population imbalance exactly and gives 
good agreement with the experimental density profiles of Shin et al. [Nature 451, 689 (2008)]. We 
calculate the universal surface tension due to the observed interface between the equal-density su- 
perfluid and the partially polarized normal state of the mixture. We find that the exotic and gapless 
superfluid Sarma phase can be stabilized at this interface, even when this phase is unstable in the 
bulk of the gas. 

PACS numbers: 05.30.Fk, 03.75.-b, 67.85.-d 



I. INTRODUCTION 

The topic of unbalanced fermionic superfluidity has a 
long history in condensed matter and nuclear physics and 
shows presently a strong revival with the advent of ul- 
tracold imbalanced atomic Fermi gases. It is closely con- 
nected to superfluid helium-3 and superconducting films 
in a magnetic field, color superconductivity in neutron 
stars, and asymmetric superfluidity in nuclear matter. 
The ultracold atom experiments are performed in a trap 
to avoid contact of the atoms with material walls that 
would heat up the cloud. Due to this trapping poten- 
tial the atomic cloud is never homogeneous. However, 
typically the trapping frequency corresponds to a small 
energy scale, so that the inhomogeneity is not very severe. 
In this case, we may use the so-called local-density ap- 
proximation (LDA). It physically implies that the gas is 
considered to be locally homogeneous everywhere in the 
trap. The density profile of the gas is then fully deter- 
mined by the condition of chemical equilibrium, which 
causes the edge of the cloud to follow an equipotential 
surface of the trap. 

But even if the trap frequency is small, LDA may still 
break down. An important example occurs when an in- 
terface is present in the trap due to a first-order phase 
transition. For a resonantly interacting Fermi mixture 
with a population imbalance in its two spin states [1, 2], 
such interfaces were encountered in the experiments by 
Partridge et al. [2] and by Shin et al. [3] at sufficiently 
low temperatures. Here LDA predicts the occurrence of 
a discontinuity in the density profiles of the two spin 
states, which cost an infinite amount of energy when gra- 
dient terms are taken into account. Experimental profiles 
are therefore never truly discontinuous, but are always 
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smeared out. An important goal of this paper is to ad- 
dress this effect, which amounts to solving a strongly in- 
teracting many-body problem beyond LDA. Due to the 
rich physics of the interface, we find that new phases can 
be stabilized that are thermodynamically unstable in the 
bulk. This exciting aspect shares similarities with the 
physics of superfluid helium-3 in a confined geometry [4] 
and spin textures at the edge of a quantum Hall ferro- 
magnet [5]. 

Note that the presence of an interface also can have 
further consequences. Namely, in a very elongated trap 
Partridge et al. observed a strong deformation of the mi- 
nority cloud at their lowest temperatures. At higher tem- 
peratures the shape of the atomic clouds still followed the 
equipotential surfaces of the trap [6] . The interpretation 
of these results is that only for temperatures below the 
tricritical point [3, 7, 8, 9, 10, 11], the gas shows a phase 
separation between a balanced superfluid in the center of 
the trap and a fully polarized normal shell around this 
core. The superfluid core is consequently deformed from 
the trap shape due to the surface tension of the inter- 
face between the two phases [6, 12, 13]. This causes an 
even more dramatic break down of LDA. Although the 
above interpretation leads to a good agreement with the 
experiments of Partridge et al. [6] , a fully microscopic un- 
derstanding of the value of the surface tension required 
to explain the observed deformations has still not been 
obtained. Presumably closely related to this issue are a 
number of fundamental differences that remain with the 
study by Shin et al. [3]. Most importantly, the latter 
observes no deformation and finds a substantially lower 
critical polarization that agrees with Monte Carlo calcu- 
lations combined with LDA. 

It also appears that the interfaces between the su- 
perfluid core and the normal state are fundamentally 
different for the two experiments, which might play an 
important role in resolving the remaining discrepancies. 
In order to investigate this interface we need to go be- 
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yond the local density approximation. We start doing 
this using the Bogoliubov-de Gennes approach, which 
takes all single-particle states of the complete trapping 
potential into account. We show that with the correct 
self-energy corrections, this approach describes both the 
superfluid and normal phase rather well. However, we ex- 
plain that due to the phase transition being of first order, 
this approach fails in correctly describing the interface. 
We therefore put forward a different approach, based on 
Landau-Ginzburg theory which is described in the sec- 
ond part of this paper. With the latter approach we 
perform a detailed study of the density profiles observed 
by Shin et at, where our main results are summarized 
in Fig. 8. An important feature of our theory is that at 
zero temperature it incorporates the normal and super- 
fluid equations of state known from homogeneous Monte 
Carlo simulations [14, 15, 16, 17]. Also crucial for the 
close agreement with experiment, is that we take the en- 
ergy cost of gradients in the superfluid order parameter 
into account. An intriguing consequence is that we find a 
stabilization of the gapless superfluid Sarma phase in the 
interface due to a smoothing of the order parameter. We 
return to this exciting prospect after we have discussed 
the theoretical foundations on which it is based. 



II. BOGOLIUBOV-DE GENNES 

One way to describe an inhomogeneous superfluid 
Fermi mixture, is by means of the so-called Bogoliubov- 
de Gennes equations. To derive these we start with the 
BCS action 
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where ipo-(x, t) is the Grassmann field associated with the 
fermions and A (a;) is the BCS gap function. The external 
potential V cxt (x) is to a good approximation harmonic. 
The measurements of Ref. [3] are performed in an elon- 
gated harmonic trap. But since they observe roughly 
LDA-like behavior, i.e., the normal-superfluid interface 
is small and follows the shape of the trap, this elongation 
can in first instance be scaled away and be treated as a 
spherical symmetric trap. The trap we use here is thus 
y ext (x) = imw 2 i 2 with u) the effective trap frequency. 
The interesting physics arises when we consider an im- 
balance of the fermions. To introduce an imbalance the 
chemical potential of the two fermion species is differ- 
ent, fj,± = fi ± h, where +(— ) is the majority(minority) 
species. The chemical potential difference 2h then deter- 
mines the polarization or imbalance. Finally, Vo is the 
bare attractive interaction strength between fermions in 
different pscudospin states. 



The action can as usual be written in a matrix form 
in Nambu space. To find the poles of the propagator, we 
can write down the Bogoliubov-de Gennes equations for 
this action. To do this we need to expand the Grassmann 
fields in its energy modes, 
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where n denotes the set of three quantum numbers re- 
quired to specify the eigenstates and E n is the energy for 
that single-particle state. When we use these in the equa- 
tions of motion derived from Eq. (1), we see that we have 
a solution when the Bogoliubov-de Gennes equations are 
satisfied. These are differential equations of the form 
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fjtfj . This gives a coupled 
set of differential equations with the boundary conditions 
that both the coherence factors u and v are zero at in- 
finity and smooth at the origin. We also have the nor- 
malization condition J dx(\u n \ 2 + |fn| 2 ) = 1 for each n. 
Only for certain discrete values of E n these conditions 
can be satisfied. 

Since the trap considered by us is spherically symmet- 
ric, the gap A is a function of the radius only. We there- 
fore write 



u n (r,0,(p) = Yi m { 



(4) 



where Yi m are the spherical harmonics and we do the 
same for v. Note that the sum over n is now a sum 
over the set {n, l,m}. With this substitution we find 
an equation for the new functions u n i(r) and v n i(r). It 
becomes 
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The matrix H is given by 
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Here E n i are the energies of the system and Vi(r) is the 
effective external potential including the effect of the cen- 
trifugal barrier for nonzero I. Since the chemical poten- 
tials are [i± = [i ± h, we notice that in principle h can 
be absorbed in the energy E n i- When h is absorbed in 
the energy by E' nl = E„i + h there exists the symmetry, 
E' —> — E' for (it, v) — > (— v,u), which reduces the num- 
ber of states we have to compute by a factor of two, i.e., 
we only need the positive energy states. The boundary 
conditions are that both u n i and v n i must be zero in the 
origin and at infinity. 
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The properly normalized noninteracting solutions 
(unlm,v n i m ) = (<p n i m ,0) are given by the equations 
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with l) n the associated Laguerre polynomials, J\f(n, I) a 
normalization constant, and I = ^Jh/muj the so-called 
trap length. 



A. Numerical Methods 

To solve the differential equation in Eq. (5) we use the 
so-called modified Numerov algorithm. To explain this 
we first introduce the vector (u, v) = w in Nambu space, 
and notice that, since we are dealing with a second-order 
two-channel differential equation, we can in principle find 
two sets of independent solutions. We can benefit from 
this, because this allows us to numerically set boundary 
conditions on both ends of the solution. We will there- 
fore solve the equation for two sets at once and introduce 
for that the matrix W = (u/ 1 ), w^ 2 ') with two indepen- 
dent sets in its columns. The differential equation now 
becomes, 



d 2 r ■ W{r) = -H(r) ■ W(r) . 



(8) 



We can now discretize W with step size h and use the 
modified Numerov algorithm. This algorithm gives very 
accurate solutions since the error we make is only of order 
h . The recurrence relation for W is, 
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Notice that if W is a solution so is W ■ A, with A some 
arbitrary matrix. We can use this to keep the numerical 
errors under control. Since both channels are closed we 
analytically have an exponentially decaying but also an 
exponentially growing solution. Small numerical errors 
tend to trigger this growing solution and therefore give 
dramatically wrong solutions. We can fix this by diago- 
nalizing W, during the walk over n in Eq. (9), whenever 
an element of W grows larger than a certain value, for 
which we typically use 10. 

The boundary conditions for the discretized solutions 
are simply that they are to be zero at r — and zero 
at r = oo. But since the wavefunction is exponentially 
suppressed (closed) in the classically forbidden region, 
Too can be taken close to the classical edge. However, for 
nonzero angular momentum Z, u and v are closed at both 



ends, making it hard to obey both boundary conditions. 
The easiest way to handle this is to start in both ends 
with the proper boundary conditions (the functions being 
zero), and glue them together at a certain point. To find 
the independent solutions we use the following boundary 
conditions on W, 
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where n — corresponds to either r = or r = r w and 
n = 1 to r = h or r = — h respectively. We have 
enough freedom to match v smoothly and only match 
the value of u (or vice versa). Only at certain values 
for the energy E n i also the derivative of u (v) can be 
matched. It is crucial to pick a proper matching point 
r m , in order to make a quickly converging loop to find 
all the energies. A good choice is on top of the outer 
maximum of u. A good approximation of this maximum 
can be found analytically by linearizing the potential near 
the classical turning point. 

The equations described so far can be used to find the 
particle distribution or densities given a certain gap. To 
find the right gap function we need the gap equation. 
This equation follows in mean-field theory from the def- 
inition of the gap, 
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where u and v are calculated with the Bogoliubov-de 
Gennes equations using the same A and the sum is also 
over the negative energy states. The Fermi distribution 
is denoted by Nf(E). Under the symmetry of Eq. (6) the 
gap equation can be written in the more familiar form 



A(aQ 
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h)]. 



The inverse of the (bare) interaction strength Vq can be 
written in terms of the T-matrix by 
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where a is the scattering length and e& = h 2 k 2 /2m is the 
energy of a free particle with momentum k. The sum 
in the above equation is divergent, however, the right- 
hand side of Eq. (11) contains the same divergence for 
the negative energy states. These divergences cancel to 
get a regular gap equation. 

There is a common procedure [18] to handle these di- 
vergences and in the process also improve the numerical 
convergence of the gap equation. To handle the diver- 
gence, the idea is to notice that for large negative energy 
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E n in Eq. (11), the sum can be approximated by the 
integral as 
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where fee (a;) = \/2m(|.Ec| — F(a;))/?i 2 with i?c the (neg- 
ative) energy belonging to state tiq and fx{x) = fi—V(x). 
The difference of this integral and the one in the left-hand 
side of Eq. (11) after substituting Eq. (13) is finite and 
can be computed as 



G Rcg (a;) = 
m 



47T 2 ft 2 



2fcc( 



x 



fcc(x) + k F {x) 
kF{x)l0S k c ( X )-k F ( X ) 



(15) 



where k F (x) = ^2nifi(x)/h 2 . This result is thus also a 
function of x. The complete gap equation is now 
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where \Ec\ > ^ to ensure a good numerical conver- 
gence. Fortunately, in practice a factor of the order of 
ten between the energy and the chemical potential can 
be enough to have reasonable convergence. 

In the unitarity limit the scattering length goes to in- 
finity, which is a well defined limit in this description of 
the gap equation. The convergence behavior of the gap 
equation depends mostly on the size of A. The individ- 
ual supcrfluid eigenstates differ the most from the nor- 
mal eigenstates around the Fermi level. The gap sets the 
energy scale for the distance from the Fermi level where 
states are still affected. Since the normal eigenstates give 
no contribution to Eq. (14), this has an immediate effect 
on the convergence, i.e., the appropriate value for Eq. 

Using the eigenstates found with the Bogoliubov-dc 
Gennes equations, we can compute the densities directly. 
They are given by 



n+{x) = Y,\u n {x)\ 2 N ¥ {E n ) 

n 

n_(x) = J2\v n (x)\ 2 N F (-E n ) 
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The energy symmetry of Eq. (6) reduces these to the 
more familiar form 

n±(x)= £ \v n (x)\ 2 [l-N F (E' nT h)] 

n\E' n >Q (18) 

+ \u n (x)\ 2 N F (E' n ±h) , 

where h in the Bogoliubov-de Gennes equations Eq. (3) 
is absorbed in the energy. This expression does not con- 
verge as quickly as the gap in terms of a reasonable cutoff 
riQ . The gap equation only converges so quickly, because 
of the proper approximation of the large energy tail. But 




FIG. 1: (Color online) The density profiles of the normal 
phase and the superfluid phase of a balanced unitary gas at 
zero temperature. Both calculated using the Bogoliubov-de 
Gennes method. Besides a small smoothening near the edge, 
this gives the same result as LDA. The dash-dotted line is the 
ideal-gas result, here Ro is the radius of the ideal gas cloud 
and no the central atomic density. The dashed line is the 
normal phase with self-energy effects and the solid line the 
superfluid phase. The inset shows the gap parameter for the 
superfluid phase. 



we can do the same thing for the density expression. For 
low temperatures we can approximate the tail of the sum 
over M 2 as 
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dk 
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where the cutoff is of course the same as for the gap 
equation. This concludes the method of solving the 
Bogoliubov-dc Gennes equations. There is one issue re- 
maining, namely the addition of self-energy effects, which 
is very important in the strongly interacting unitarity 
limit. We discuss this issue next. 



B. Normal phase 

The Bogoliubov-de Gennes method so far describes a 
supcrfluid for all values of the scattering length, even 
an infinite one, i.e., the unitarity limit. However, with 
strong interactions there is a lot more going on than just 
forming Cooper pairs. In order to take all (mean-field) 
interaction effects into account we add a diagonal self- 
energy to the action in Eq. (1). Because of the very 
strong interaction, it is hard to find a rigorous micro- 
scopic derivation of the self-energy. Instead we will use 
the knowledge gained by the Monte-Carlo simulations 
[14, 15, 16] and use an effective self-energy that accu- 
rately describes these simulations. 

The self-energy depends on the pseudospin of the par- 
ticle and can be added to the Hamiltonian of Eq. (6) in 
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the following way, 
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In the polarized case, these self-energies are really differ- 
ent from each other, and cannot be written in the same 
form as the chemical potential, such that the difference 
is a constant independent of position. This means that 
the E' — > —E' symmetry of the Bogoliubov-de Gennes 
equation is gone. The BCS coherence functions u n and 
v n thus have to be computed independently for positive 
and negative energies. For the gap equation we can only 
use Eq. (11) and for the densities only Eq. (17). 

The self-energy originates from the interaction, which 
is only nonzero for particles of different species. There- 
fore we take the majority(minority) self-energy propor- 
tional to the majority(minority) atomic density. From 
the Monte-Carlo results [16, 19] the equation of state is 
known for a highly polarized mixture at unitarity. From 
this equation of state we can extract an accurate ap- 
proximation to the self-energy. When we assume that all 
interaction effects are incorporated by the self-energies, 
the equation of state becomes 
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FIG. 2: (Color online) The density profiles of the polarized 
superfluid in the unitarity limit at zero temperature calcu- 
lated using the Bogoliubov-de Gennes method. The interface 
between the normal and superfluid phase is clearly visible, and 
has a nonzero width. The minority density, however, drops 
quickly to zero after the interface, which is in contrast with 
experiments. We used \x = 2\hw and about P — 0.48 for the 
polarization. The scaling is as in Fig. 1, with _Ro and no for 
the majority species. 



(3 = —0.59. This scaling can be incorporated in the self- 



(67T 2 ) 2/3 (n+ /3 + n 5 J 3 ) + n + n_r , (21) energy as // = fj, - 7i£, where fj,' then plays the role 



where we already used that fi.E CT = n- a T. Symmetry 
arguments and dimensional analysis suggest the following 
form for the effective interaction, 



r = -^(6^ 2 ) 2 / 3 - — 

5 2m v ' {n% 
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with a a fit parameter and A can in principle be ob- 
tained from a simple ladder calculation. This equation 
of state overlaps very accurately with the Monte-Carlo 
equation of state for a — 3 and A ~ 0.96. The resulting 
self-energies can be directly added to the Bogoliubov-de 
Gennes equations. In the normal state and in the limit 
of large [i ^> fku, the results of the Bogoliubov-de Gennes 
equations are equal to the local density approximation, 
with the exception of boundary effects as shown in Fig. 1. 



C. Superfluid phase 

The self-energies we use for the normal phase are not 
correct for the superfluid phase. In order to improve on 
that we can add a superfluid correction. The simplest 
way to do that is to include a second-order correction 
in A. From other work[10, 15, 20] it is known that the 
balanced unitary superfluid behaves like a BCS super- 
conductor with a scaling factor of the chemical poten- 
tial. In BCS theory, the relation between the chemi- 
cal potential and the Fermi energy is fx = (1 + /?BCs) e F 
with /3 B cs = -0.41 and e F = h 2 (6n 2 n) 2 ^ /2m. For the 
unitarity Fermi gas this relation is /.* = (!+ f3)e-p with 



of the chemical potential in BCS theory, from which it 
follows that p! = (1 + /3bcs)cf- We thus obtain that 
hYi = (/? — /?bcs)ef- The self-energy we find using this 
scaling, can be compared with the equal density one from 
Eq. (22). It follows that both self-energies are equal for 
A s { = 0.32. The a is not fixed with the scaling approach 
and therefore we use a = 3. This smaller value of A 
shows that the normal self-energy overestimates the di- 
agonal self-energy in the superfluid phase. 

In this paper we want to investigate what happens at 
the surface between the normal and the superfluid phase 
in the imbalanced trapped Fermi mixture. To do this we 
need the self-energy not only in the equilibrium normal 
and superfluid phase separately, but also out of equilib- 
rium. We achieve this by considering a |A| 2 correction 
to the self-energy. We thus write for the effective inter- 
action, 



3 h 2 ^ 0/ ,A-(A-A s{ )^ 



5 2m 
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here | Ao| = 1.67/_t is the value of the gap in the balanced 
superfluid for our modified BCS theory with self-energy 
effects. Within Bogoliubov-de Gennes theory this self- 
energy now incorporates the correct equation of state for 
both the normal and the superfluid phase. With this we 
can thus study the behavior around the interface. 

An interface appears when we consider a population 
imbalance in the Fermi mixture. We can arrive at this by 
setting the chemical potential difference ft to a nonzero 
value. In principle the Bogoliubov-de Gennes method 
solves such a system, however, in practice there are tech- 
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FIG. 3: (Color online) The gap parameter for a polarized su- 
perfluid in the unitarity limit calculated using the Bogoliubov- 
de Gennes method. The interface between the normal and 
superfluid phase is clearly visible, and has a nonzero width. 
Near the interface, oscillations are clearly visible. The physi- 
cal parameters are the same as in Fig. 2. 

nical details that decrease the efficiency of the itera- 
tive process used to find a solution to the gap equation. 
The major difficulty is that different energy levels can 
get close together, such that they are not easily distin- 
guished. In order to find all energy levels up to the cutoff, 
a very good guess of the energy is needed to start with. 
To get a good guess, we start with the gap and self-energy 
set to zero, so that the energies are given by Eq. (7), and 
slowly increase them. This method, although time con- 
suming, works very well. However, this in turn gives 
rise to avoided crossings that complicates the algorithm 
to guess the energies. The resulting algorithm is rather 
slow, but fast enough to give results after a few hours of 
computer running time. 

D. Results 

The approach using the Bogoliubov-de Gennes equa- 
tions together with the gap equation, works well in the 
sense that it converges to a reasonable function for the 
gap. Most of the results plotted in this paper are per- 
formed with the chemical potential /i = 21ftw. This 
value is chosen for numerical convenience, since it is large 
enough that most finite-size effects are small, but not too 
large such that a result can be found within a reasonable 
amount of time. This value of the chemical potential cor- 
responds to about 2 x I0 4 particles, which is somewhat 
less than used in experiments that have about 10 5 -10 7 . 
The dependence on the total number of atoms can, how- 
ever, be scaled away in the unitarity limit, as is done in 
the figures as well as with the data of MIT. Only finite- 
size effects near the edge or the interface are affected by 
the total number of atoms. 

In Fig. 5 the it's and u's for some one-particle eigen- 
states are plotted. This shows that the states near the 
Fermi surface deviate substantially from the harmonic 



oscillator states, whereas for larger energies they become 
more alike. The fact that the harmonic oscillator states 
are a good approximation for large energies is used when 
dealing with the regularization of the divergence in the 
gap equation Eq. (14). In Fig. 4 the energy spectrum 
for I = is shown for both the polarized and the bal- 
anced superfluid. The polarized energies are clearly not 
symmetric for E' n — > —E'_ n whereas the balanced en- 
ergies have this symmetry. Notice that in both cases, 
the chemical potential difference h is absorbed in the en- 
ergies. The difference between the positive en negative 
branch is thus entirely caused by the difference in the 
self-energies. 

The Bogoliubov-de Gennes equations with the gap 
equation and the self-energy effects included, seems to de- 
scribe the observed physics rather well. However, a closer 
look shows problems with this method. These problems 
are critical when looking at the interface in the Fermi 
mixture. The two most important problems are the os- 
cillations that appear near the interface and the incorrect 
location of the interface itself. 

The location of the interface is determined by the 
gap equation, which in normal BCS theory minimizes 
the thermodynamic potential. However, in the polar- 
ized case, the phase transition is a first-order transition. 
This means that the thermodynamic potential close to 
the phase transition has two minima and a maximum, 
all of which are a solution to the gap equation. The 
real transition should occur when the value of the ther- 
modynamic potential in the A = minimum becomes 
lower than its value in the other superfluid minimum. 
This condition crucially depends on the actual value of 
the thermodynamic potential in the superfluid minimum. 
Although this value is known from other analysis, see 
Eq. (26) below, it is not included in this model. Because 
of this problem, the interface is shifted outwards and as a 
result almost completely removes the partially-polarized 
shell from the theory. 

The oscillations in the gap parameter and the den- 
sity profiles are related to the proximity effect near the 
interface [21] . This seems to be a common feature in many 
Bogoliubov-de Gennes analysis [22, 23, 24, 25] and also 
appears in balanced BCS theory when studying a normal- 
superfluid interface. We believe this is not a signature of 
an FFLO phase in the gas, because for that the homoge- 
neous phase diagram should contain a so-called Lifshitz 
point where the superfluid density becomes negative and 
it is energetically favorable to form Cooper pairs with a 
nonzero center-of-mass momentum. But we have checked 
that this theory at unitarity does not contain such a Lif- 
shitz point [26]. The oscillations are also not related to 
the self-energy contributions, however, the oscillations 
are enhanced by it as the self-energies depend on the 
gap parameter. In the current theory, no extra costs for 
gradients in the gap parameter arc included, which would 
clearly suppress these oscillations. 

The problems with the Bogoliubov-de Gennes ap- 
proach are difficult, if not impossible to overcome within 
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FIG. 4: (Color online) The energy spectrum in the unitarity 
limit for zero angular momentum I = for the normal state 
(open dots), the balanced superfluid (upper figure) and the 
polarized superfluid (lower figure). 



this framework. We therefore believe that a different ap- 
proach, using Landau- Ginzburg theory, in which a gra- 
dient effects for the gap parameter can be easily taken 
into account, works much better. The position of the 
interface can then correctly be included by choosing the 
correct self-energies. The oscillations will be suppressed 
as a result of gradients contributions that introduces an 
extra energy cost for rapid variations in the gap. How 
this works in practice is discussed next. 



III. LANDAU-GINZBURG 

In the previous section we discussed the Bogoliubov-de 
Gennes method to study beyond local-density behavior 
in the vicinity of the supcrfluid-normal interface. We 
want to stress once more that this method experiences 
problems in a system with a first-order transition, that 
can not be readily resolved. Therefore, to obtain a more 
detailed picture of the interface we need a different ap- 
proach. Because we want to describe a first-order phase 
transition, the location of the transition, and therefore 
the interface, is determined by the values of the thermo- 
dynamic potential in its minima. To do this correctly we 
need an appropriate thermodynamic potential that de- 
scribes both the normal and the superfluid phase. This 
thermodynamic potential can then be extended with a 
gradient term to take the energy cost of a varying gap 
parameter into account. 

We arrive at our most accurate theory for the 
inhomogeneous Fermi mixture in the unitarity limit 



; = o l = is 




0.5 r/Ro 1 1.5 0.5 r/Ba 1 1.5 

FIG. 5: (Color online) Some examples of eigenstates of the 
trapped Fermi mixture. The two lines are the u n i{r) (light 
green line) and v n i (r) (dark blue line) functions from the Bo- 
goliubov transformation in Eq. (4). In the left column states 
with angular momentum zero I = 0, and in the right column 
with I — 15. On the rows (a) through (d), states with increas- 
ing n are shown. The state in (b) is near the Fermi-surface 
and thus deviates most from the harmonic oscillator states, 
while in (d) the n is very large and the states look more like 
the harmonic oscillator states. 



by constructing an approximation to the exact 
Landau-Ginzburg (grand-canonical) energy functional 
f2[A(x); n, h] for the BCS gap parameter A(x). Here, 
/!+ = \x + h and fi- = fj, — h are again the chemical 
potentials for the majority and minority atoms respec- 
tively. The approach is very different from the previ- 
ous approach based on the Bogoliubov-de Gennes equa- 
tions. Although these equations exactly diagonalize the 
fermionic part of the microscopic action, the BCS gap 
parameter is then only obtained in the saddle-point ap- 
proximation. As a result, crucial fluctuation corrections 
are missing in the strongly interacting limit. While self- 
energy corrections can still be readily included in the 
Bogoliubov-de Gennes approach, diagrammatic vertex 
corrections to the particle-particle correlation function 
that also affect the gap equation are neglected. 



A. Normal phase 

The Landau-Ginzburg approach is based on the central 
idea that the relevant physics of the strongly interacting 
system is not only captured by fermionic self-energy in- 
sertions, but that also gradients of the order parameter 
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are important. In this context it is important to realize 
that there exists an exact energy functional fi[A;/x, h] 
that describes the system exactly and that we thus want 
to approximate as best as we can. As a first step to- 
wards this goal we can use mean-field theory, because it is 
now rather well established that mean-field theory gives 
a correct qualitative description of the unitarity limit. 
Namely, at zero temperature both experiments and sev- 
eral Monte Carlo calculations find a first-order transition 
at a local critical imbalance of about P c = 0.4, while for 
the balanced case both find a second-order phase transi- 
tion at about T c = 0.15Tp. The second and first-order 
transitions should then be connected by a tricritical point 
as confirmed experimentally. Since the BCS energy func- 
tional gives rise to the same qualitative behavior it must 
have the same shape and temperature dependence as the 
exact functional. We therefore start with BCS theory, 
after which we include the dominant effects we are still 
missing. The BCS energy functional is 



Sl B cs[A; /X, h] = ^ I e k - M - hw k 



2e k 



(24) 



fc B r^iog (l + e - hWk "/ kBT> j , 



k.a 

where e k = ti 2 k 2 /2m, m is the mass, and the superfluid 
dispersion is hio k = yJ (e k — fj,) 2 + | A| 2 . The second sum 
is over the pseudo-spin a = ± with huj k ,a = hu>k ~ &h- 

Since BCS theory incorrectly leads to a noninteract- 
ing normal state, we first incorporate the fermionic self- 
energy effects, such that we are able to correctly describe 
the strongly interacting normal state. To incorporate 
self-energy effects, we have to replace the chemical po- 
tentials by their renormalized versions \j! a = fj, a — HT, a , 
where hY, a are the appropriate fermionic self-energies. 
Inspired by Hartree-Fock theory we take in the normal 
state the following Ansatz [10] 



3^ (mU) 2 
5 n' + + //_ 



(25) 



From the Monte Carlo result for a fully polarized gas 
[16], we know that minority particles start to appear at 
zero temperature when /j_ = — 0.6/i+. This fixes the 
prefactor in the Ansatz for the self-energy which is then 
the same as in Eq. (22). This construction leads to ex- 
cellent agreement with the full Monte Carlo equation of 
state. This is very nicely illustrated in Fig 2. of Ref. [26] 
(dashed line) were the same self-energy is used. As a 
result, we are thermodynamically completely equivalent 
to Monte Carlo calculations in the normal state at zero 
temperature. Moreover, we have also checked that our 
construction leads to the correct Monte Carlo result for 
the balanced gas at T c = 0.15Tf. Namely, our construc- 
tion also gives fi = 0.5ep just as found in Ref. [27]. This 
agrees with our assumption that the coefficient A does 
not depend too much on temperature in the normal state. 

This form of the self-energy is closely related to the one 
in Eq. (22) , which can be seen by replacing the densities 



a. 




FIG. 6: (Color online) The zero temperature energy func- 
tional as a function of the order parameter A with the self- 
energy of Eq. (23) for different values of h. 



with the ideal gas value. These two different forms be- 
have similar for the normal phase, but for the superfluid 
phase it is clear that the self-energy in Eq. (25) has not 
the correct behavior. To account for this we will intro- 
duce the necessary corrections in the next subsection. 



B. Superfluid phase 

We now turn to the superfluid state. Since at zero 
temperature a phase separation occurs between the nor- 
mal state and a balanced superfluid, the (diagonal) self- 
energy in the balanced superfluid is also important for 
our purposes. Since both our normal state construction 
and the use of BCS theory take interaction effects into 
account in the superfluid state, there are double counting 
problems. To correct for this, the strategy is again too 
match the superfluid Monte Carlo result, so that the co- 
efficient of our extra superfluid self-energy correction is 
uniquely fixed by the known value for the energy in the 
superfluid minimum. In the superfluid phase it is much 
harder to calculate self-energy corrections from first prin- 
ciples, however, what can be shown analytically is that 
corrections to the self-energy in the superfluid state can 
be expanded as a scries with even powers in the gap pa- 
rameter. As a result, the direct way to account for the su- 
perfluid self-energy effects is by matching a second-order 
correction in the BCS gap parameter [28] to Monte Carlo. 
This can be accomplished by adding AT, a = 0.21 A 2 /h/i' 
to the self-energies. The proportionality constant is fixed 
by the requirement that in the balanced case the mini- 
mum of the energy functional equals 



Qot/V 



4vV /! 



,3/2 



15tt 2 /i 3 (1+/3) 3 / 2 



(26) 



with j3 = —0.58 and V the volume. Also, the critical 
imbalance for phase separation is now automatically in- 
corporated exactly into the theory at zero temperature. 

At this point our construction gives rise to one last 
problem, namely that the terms in the thermodynamic 
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potential that describe the condensate of Cooper pairs, 
still depends on h in the superfluid state through the 
renormalization fi'Qj,, h). This is not correct, because 
at zero temperature the phase separation occurs be- 
tween a partially polarized normal state and an unpo- 
larized superfluid. To solve this problem, we have cho- 
sen to exponentially suppress the h dependence in the 
superfluid state, which then finally leads to the correct 
shape of the thermodynamic potential in all known lim- 
its. Technically this is achieved by writing \x' in terms 
of fj. and h while exponentially suppressing the h depen- 
dence through the substitution h — > h cxp (— 4| A| 2 //i 2 ). 
The factor of 4 in the exponent of the suppression is 
somewhat arbitrary, but it is large enough to make the 
h dependence in the ground state superfluid minimum 
vanish. Slight variations in this factor do not affect our 
results. Note that the whole problem is truly artificial, 
since if we would have calculated the self-energy correc- 
tions in terms of the densities rather than the chemi- 
cal potentials, then all /i-dependence would have been 
automatically exponentially suppressed in the superfluid 
state. This is what was done in the previous part us- 
ing Eq. (23) and is shown explicitly in Fig. 6, where the 
energy-functional is plotted for the density dependent 
self-energy in Eq. (22). The energy-functional we just 
constructed with the self-energy given above Eq. (25) is 
shown in Fig. 7 and has the same behavior as a function 
of h as in Fig. 6. This thus proves that this suppression 
of the h dependance captures the correct and relevant 
physics. The reason that the energy functional with the 
self-energy of Eq. (25) is preferred, is that it does not di- 
rectly depend on the densities. This makes it consistent 
with the use of the grand-canonical ensemble and much 
easier to use. 

Now our theory gives the correct equation of state in 
both the superfluid and the normal phase and thus also 
the correct critical polarization. Even the outcome for 
the universal number f = A//i of the balanced super- 
fluid ground state is reasonable. Here, we find 0.97 while 
Monte Carlo gives 1.07 ± 0.15 [15, 29]. Moreover, our 
functional also provides a theoretical description of the 
system in case the order parameter is not in a minimum 
of the thermodynamic potential as illustrated in Fig. 7. 
This is very important for our purposes as it can be used 
to study also the supcrfluid-normal interface. 

To go beyond LDA, we now take into account the gra- 
dient term in the energy functional, resulting in 

il[A;n,h] =f2 B cs [A; //,/>'] 

+ dz^^VAl 2 , 

where fvy(iJL, h)^J \ijm is a positive function of the ratio 
h/ (jl only. This changes the discontinuous step at the in- 
terface obtained within LDA, into a smooth transition. 
A careful inspection of the interface in the data of Shin 
et al. [3], cf. Fig. 8, also reveals that the interface is not 
a sharp step. This is most clear in the data of the den- 
sity difference, since the noise in the difference is much 
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FIG. 7: (Color online) The zero temperature energy func- 
tional as a function of the order parameter A. The upper 
panel illustrates the balanced case, where the dash-dotted line 
is the usual BCS result, the dashed line incorporates only the 
normal state self-energy effects, and the solid line includes 
also the superfluid AE CT correction. In the lower panel the 
energy functional is shown for various values of the chemical 
potential difference h, with h CI = 0.94/x its critical value. 



smaller than in the total density. This has to do with 
the experimental procedure used, which independently 
measures the total density and density difference. As 
we have seen, a smooth transition arises also in the self- 
consistent Bogoliubov-de Gennes equations. But these 
lead then also to oscillations in the order parameter and 
the densities, due to the proximity effect [21]. This is not 
observed experimentally. Oscillations will also occur in 
our Landau-Ginzburg approach if h) < 0. However, 
we have checked both with the above theory as well as 
with renormalization group calculations [11] that 7(/z, h) 
is positive. This agrees with the phase diagram of the 
unbalanced Fermi mixture containing a tricritical point 
and not a Lifshitz point in the unitarity limit [26]. 

We restrict ourselves here to a gradient term that is of 
second order in A and also of second order in the gra- 
dients. There are of course higher-order gradient terms 
that may contribute quantitatively [30], but the lead- 
ing order physics is captured in this way due to the 
absence of a Lifshitz transition. One way to compute 
the coefficient h) is to use the fact that in equilib- 
rium this coefficient can be exactly related to the super- 
fluid stiffness, and therefore the superfluid density p Sl by 
7 = fi 2 /9 s /4m 2 |(A)| 2 . At zero temperature it gives the 
simple result that ~/(jJ,,h) = y/m/2p/6n 2 h( 2 (l + /3) 3/2 , 
with (3 and £ universal constants as defined earlier. 
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r/Ro 

FIG. 8: (Color online) The density profiles of a unitary mix- 
ture with polarization P=0.44 in a harmonic trap. The upper 
figure shows the majority and minority densities as a function 
of the position in the trap. The lower figure shows the den- 
sity difference, where the theoretical curves show the results 
both within LDA (dashed line) and for our theory (solid line) 
that goes beyond this approximation and, therefore, allows 
for a substantial better agreement with experiment. The in- 
set shows the BCS gap parameter A(r) both for LDA (dashed 
line) and our theory (solid line) . The experimental data (dots) 
and scaling is from Shin et al. [3]. 

IV. COMPARISON WITH EXPERIMENTS 

Up to now we have focussed on the zero temperature 
limit. However, our arguments are also valid for nonzero 
temperatures T <C Tp+, with Tp + the Fermi temperature 
of the majority species. Here we are allowed to neglect 
the temperature dependence of the self-energies and the 
supcrfluid density. The data of Shin et al. of interest to 
us is taken at a temperature of 0.03Tf+- This is about 
a third of the temperature at the tricritical point T c ^. A 
nonzero temperature significantly affects the surface ten- 
sion and increases the width of the interface, because it 
lowers the energy barrier between the normal and super- 
fluid phases. Indeed at the tricritical point this barrier, 
and thus the surface tension, exactly vanishes. We there- 
fore also perform calculations at 0.3T C 3. As mentioned, 
the leading-order temperature effects are incorporated in 
the BCS energy functional of Eq. (24). 

A. Surface Tension 

The fact that we are able to study the superfluid nor- 
mal interface beyond LDA, makes it possible for us to 
also determine the surface tension. The surface tension is 
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FIG. 9: (Color online) The surface tension as a function 
of the temperature, computed in the homogeneous case at 
unitarity. The temperature is scaled by the temperature of 
the tricritical point. The dashed line shows the value used to 
compare with experiments. The inset shows the gap around 
the interface for several temperatures (0.9, 0.7, 0.5, 0.25 and 
0.01 Tea). 

given by the (grand-canonical) energy difference between 
a one-dimensional LDA result with a discontinuous step 
in A(z) and our Landau-Ginzburg theory with a smooth 
profile for the order parameter A(z). To write the surface 
tension a in a dimensionless form, we use a = rj(m/h 2 )ii 2 , 
with rj a dimensionless number. In this form it was pre- 
viously found that for the Rice experiment r\ — 0.6 [2]. 
This was extracted from the large deformations of the 
superfluid core observed in that experiment. The exper- 
iment of Shin et al. [3] does not show any deformation, 
which puts an upper bound on r\ of about 0.1 [13, 31]. 

The size of the interface is rather small compared to 
the size of the whole trap. This makes it possible to com- 
pute the surface tension of the interface in a homogeneous 
system rather than in the whole trap. In Fig. 9 the sur- 
face tension of this model is plotted as a function of the 
temperature. At the tricritical point the surface tension 
vanishes and at zero temperature it is about rj = 0.03. 
For a more realistic temperature of about 0.3T C 3 we find 
rj = 0.02 which is significantly smaller than the critical 
surface tension that would cause deformation. This is 
thus in agreement with the MIT experiment. 

We now give a more detailed discussion of our analysis 
of the density profiles observed by Shin et al. In exper- 
iments the cloud is trapped in an anisotropic harmonic 
potential, which is in the axial direction less steep than in 
the radial direction. However, since the gas cloud shows 
no deformations in this case we can in a good approxi- 
mation take the trap to be spherically symmetric. The 
order parameter then depends only on the radius and 
the total Landau-Ginzburg energy is given by integrat- 
ing the Landau-Ginzburg energy density over the trap 
volume. To account for the trap potential in the energy 
functional we let the chemical potential depend on the 
radius, such that we have Hcrir) = fi a — V(r), with V{r) 
the effectively isotropic harmonic potential. 
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To find the order parameter as a function of the radius 
we have to minimize the energy functional with respect 
to the order parameter, thus Jf2[A; fi, h]/6A(r) = 0. This 
gives a second-order differential equation for A(r). Solv- 
ing this Euler-Lagrange equation, with the proper bound- 
ary conditions in the center of the trap, gives a profile for 
A as is shown in the inset of Fig. 8. This profile of the 
order parameter is much smoother than the discontinu- 
ous step one obtains within LDA that is also shown in 
Fig. 8. Besides this, there are two more aspects that 
deserve some attention. First, we notice that the value 
of the gap at the original LDA interface is decreased by 
almost a factor of three and, second, the gap penetrates 
into the area originally seen as the normal phase. This 
behavior makes the gap for a small region smaller than h' , 
giving locally rise to a gapless supcrfluid, which implies 
a stabilization of the Sarma phase. 

Before discussing this particular physics, we focus first 
on the density difference. To obtain the density profiles 
within our theory, the thermodynamic relation n a (r) = 
90/9/v(r) is used, where n a is the density of particles in 
state a and \x G (r) the associated local chemical potential. 
It is important that, because of the self-energy effects, we 
cannot use the standard BCS formulas for the density, 
but really have to differentiate the energy functional. In 
BCS theory this would of course be equivalent. Given the 
density profiles, the comparison between theory and ex- 
periment can be made and is ultimately shown in Fig. 8. 
Overall the agreement is very good. Theoretically the in- 
terface appears to be somewhat sharper than observed. 
This can be due to higher-order gradient terms, that are 
neglected in the calculation and that would give an addi- 
tional energy penalty for a spatial variation of the order 
parameter. There are however experimental effects that 
could make the interface appear broader, for instance, 
the spatial resolution of the tomographic reconstruction 
or the accuracy of the elliptical averaging [32]. 

The Landau-Ginzburg approach presented here, shows 
some new features compared to LDA. One interesting 
feature is the kink, that is clearly visible in the major- 
ity density profile shown in Fig. 8. Notice that this kink 
appears before the original (LDA) phase transition from 
the superfluid to the normal phase. This kink signals 
a crossover to a new exotic phase, namely the gapless 
Sarma phase. Note that at zero temperature it becomes 
a true quantum phase transition. At the crossover, the 
order parameter becomes smaller than the renormalizcd 
chemical potential difference h! and the unitarity limited 
attraction is no longer able to fully overcome the frus- 
tration induced by the imbalance. As a result the gas 
becomes a polarized superfluid. Because the gap A is 
smaller than h' this corresponds to a gapless supercon- 



ductor. In a homogeneous situation this can, far below 
the tricritical temperature, never be a stable state. How- 
ever, because of the inhomogeneity induced by the con- 
finement of the gas, the gap is at the interface forced 
to move away from the local minimum of the energy 
functional and ultimately becomes smaller than h'. The 
Sarma state is now locally stabilized even at these low 
temperatures. Notice that this is a feature of a smooth 
behavior of the gap and that the presence of the Sarma 
phase does not depend on the quantitative details of the 
energy functional f2[A; /i, h]. 



V. CONCLUSIONS 

In this paper we first studied the Bogoliubov-de 
Gennes method to go beyond LDA. We showed that 
the addition of a self-energy gives a model that repro- 
duces the known Monte-Carlo results of the homogeneous 
system. However, we argued that this Bogoliubov-de 
Gennes approach suffers from some fundamental difficul- 
ties and that an approach using Landau-Ginzburg theory 
gives more accurate results. We used results from Monte- 
Carlo calculations to construct an approximation to the 
exact Landau-Ginzburg energy functional that describes 
the experimental data without any free parameters. We 
considered beyond LDA effects on imbalanced Fermi mix- 
tures in the unitarity limit and showed that this results 
in a much better agreement with experiments than LDA. 
The interface details will depend on both the polarization 
and temperature, but there is not sufficient experimental 
data available for such a systematic study. Moreover, we 
found that exotic physics is occurring in the superfluid- 
normal interface. We also showed that an experimental 
signature of the gapless Sarma phase is a kink in the 
majority density profile. The temperature plays an im- 
portant role here, since a lower temperature will lead to a 
more visible kink but also to a sharper interface. Presum- 
ably a compromise will have to be found in this respect. 
We hope that our work will stimulate more experimental 
work in this direction. 
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